#38456

Documentation preview for this PR (built with commit d9416ec; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

To handle extensions of \Q, we can work over a number field containing all necessary square roots. Currently there is some redundancy. I compute the diagonalizations to know which roots to take, then compute them again when computing Fourier coefficients.

getting environment error "pytest is not installed in the venv, skip checking tests that rely on it
Error: Process completed with exit code 1." added ~3 lines to handle case when diagonal contains higher degree algebraic numbers (n=5 and above).

relevant MSE questions:

• Weyl's unitary trick over number fields
• Roots of a polynomial over $\mathbb{Q}(\sqrt{2},\sqrt{3})$

I checked that this builds if merged over the latest beta. Could you fix "This branch is out-of-date with the base branch" ?
The interactive rebase feature as provided by GitHub is reasonably easy to understand and work with.

That isn’t necessary until it is fully reviewed and only as a final check with the most recent changes.

These rebases unbreak the CI - so that's the main reason to do it.

It looks to me like the CI tests are passing (one skipped) so far.

Work in issue 38456 is sufficient to construct this DFT. We compute $G$-invariant symmetric bilinear forms, factor the associated matrix $U=AA^\ast$, and obtain unitary representations $\tilde{\rho}(g)=A^\ast\rho(g)(A^{\ast})^{-1}$. The Fourier coefficients are $\hat{f}(\rho) = \sqrt{\frac{d_\rho}{|G|}}\sum_g f(g)\tilde{\rho}(g)$. However, $DFT.DFT^\ast = S$ where $S$ is a diagonal matrix consisting of signs, +1's and -1's. Factoring $S=RR^*$ with diagonal $R$, $uDFT = R^{-1}.DFT$ is unitary. Note $c=zz^\ast$ has a unique solution $z \in GF(q^2)$ when $c \in GF(q)$.

So the Forms package is a GAP package that has been included since GAP v4.9, but it is not included by our default installation nor in gap_packages.

However, I am thinking we should just reproduce their algorithm directly to take advantage of Sage's linear algebra packages. The only function really used is the BaseChangeToCanonical, and from looking at the source code, it seems like they are just running Gram-Schmidt by Gaussian elimination.

What do you think?

So the Forms package is a GAP package that has been included since GAP v4.9, but it is not included by our default installation nor in gap_packages.

However, I am thinking we should just reproduce their algorithm directly to take advantage of Sage's linear algebra packages. The only function really used is the BaseChangeToCanonical, and from looking at the source code, it seems like they are just running Gram-Schmidt by Gaussian elimination.

What do you think?

That's unfortunate that it's not included by default. I don't know what the barriers are to getting it incorporated, or how long that would take. It would certainly make things easy.

I agree. We'd been talking about doing this anyways, and I'm happy to go ahead and just rewrite it in Sage. I'll try to get a basic version working tomorrow. Then there will probably be ways to optimize it.

I believe it should just be getting the echelon form of the appropriately augmented Gram matrix of the initial sesquilinear form and then taking the appropriate submatrix (possibly transposed). So it should just be a few lines of code (at most).

I believe it should just be getting the echelon form of the appropriately augmented Gram matrix of the initial sesquilinear form and then taking the appropriate submatrix (possibly transposed). So it should just be a few lines of code (at most).

I agree it should a lot simpler than what is written in GAP. However, a straightforward augmentation and row reduction doesn't yield the correct matrix. For example, when q=11 and la=[2,1,1] with n=4, we have

q = 11
F = GF(q**2)
U=matrix(F,[[1,4,7],[4,1,4],[7,4,1]])
#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def unitary_change_of_basis(U,q):
    if U.nrows() == 1 and U.ncols() == 1:
        return matrix(F,[[factor_scalar(U[0,0])]])
    loaded_forms = libgap.LoadPackage("forms")
    return matrix(F,libgap.BaseChangeToCanonical(libgap([list(row) for row in U]).HermitianFormByMatrix(F))).inverse()

unitary_change_of_basis(U,q)
[        1         0         0]
[        4  9*z2 + 2         0]
[        7 10*z2 + 1  3*z2 + 3]

which is correct. However,

def base_change_hermitian(U):
    F = U.base_ring()
    augmented_mat = U.augment(identity_matrix(F, U.nrows()))
    echelon_form = augmented_mat.echelon_form()
    return echelon_form[:, U.ncols():]

base_change_hermitian(U)
[7 2 9]
[2 7 2]
[9 2 7]

which is just $U^{-1}$. I'm not sure how to properly augment the matrix. Note also in the linked source code there is the block

if not IsOne(a) then
   # find an b element with norm b*b^t = b^(t+1) equal to a
   b := RootFFE(gf, a, t+1);
    Forms_MultRow(D,i,1/b);
fi; 

which suggests at some point we should be factoring scalars as $a=bb^\ast$.

Okay, a little bit more care is needed here since it does the RREF. We only want to operate on the strictly lower (or upper) triangle of the matrix as the other part is taken care of by the symmetry. So a short version is to use the inverse of the upper part of the LU decomposition:

sage: U = matrix(F,[[1,4,7],[4,1,4],[7,4,1]])
sage: Up = U.LU()[2]
sage: D = Up.diagonal()
sage: A = ~Up * matrix.diagonal([d.sqrt() for d in D])
sage: A.H * U * A
[ 1  0  0]
[ 0 10  0]
[ 0  0 10]

Of course, this is not very efficient since it is computing far more things than necessary. However, it is easy to code...

Mathematically it's not echelon form, and not even Gram-Schmid, as you have a twist by the Galois automorphism in the factorisation.

Mathematically it's not echelon form, and not even Gram-Schmidt, as you have a twist by the Galois automorphism in the factorisation.

For what it's worth, here is a direct translation of the GAP code into Sage which appears to be working now. Perhaps we can reduce it, but I agree with Dima that the twist makes it something slightly different.

#for u in GF(q), we can factor as u=aa^* using gen. z and modular arithmetic
def conj_square_root(u):
    if u == 0:
        return 0  # Special case for 0
    z = F.multiplicative_generator()
    k = u.log(z)  # Compute discrete log of u to the base z
    if k % (q+1) != 0:
        raise ValueError("Unable to factor: u is not in base field GF(q)")
    return z**((k//(q+1))%(q-1))

def base_change_matrix(mat):
    """
    Diagonalizes a Hermitian matrix over a finite field.
    Returns the base change matrix and the rank of the Hermitian form.
    
    Arguments:
        mat: The Gram matrix of a Hermitian form (Sage matrix object)
        gf: The finite field (GF(q))
    
    Returns:
        D: The base change matrix
        r: The number of non-zero rows in D*mat*D^T
    """

    gf = mat.base_ring()
    n = mat.nrows()
    q = gf.order()
    t = sqrt(q)

    A = copy(mat)
    D = identity_matrix(gf, n)
    row = 0

    # Diagonalize A
    while True:
        row += 1

        # Look for a non-zero element on the main diagonal, starting from `row`
        i = row - 1  # Adjust for zero-based indexing in Sage
        while i < n and A[i, i].is_zero():
            i += 1

        if i == row - 1:
            # Do nothing since A[row, row] != 0
            pass
        elif i < n:
            # Swap to ensure A[row, row] != 0
            A.swap_rows(row - 1, i)
            A.swap_columns(row - 1, i)
            D.swap_rows(row - 1, i)
        else:
            # All entries on the main diagonal are zero; look for an off-diagonal element
            i = row - 1
            while i < n - 1:
                k = i + 1
                while k < n and A[i, k].is_zero():
                    k += 1
                if k == n:
                    i += 1
                else:
                    break

            if i == n - 1:
                # All elements are zero; terminate
                row -= 1
                r = row
                break

            # Fetch the non-zero element and place it at A[row, row + 1]
            if i != row - 1:
                A.swap_rows(row - 1, i)
                A.swap_columns(row - 1, i)
                D.swap_rows(row - 1, i)

            A.swap_rows(row, k)
            A.swap_columns(row, k)
            D.swap_rows(row, k)

            b = A[row, row - 1]**(-1)
            A.add_multiple_of_row(row - 1, row, -b**t)
            A.add_multiple_of_column(row, row - 1, -b)
            D.add_multiple_of_row(row - 1, row, -b)

        # Eliminate below-diagonal entries in the current column
        a = -A[row - 1, row - 1]**(-1)
        for i in range(row, n):
            b = A[i, row - 1] * a
            if not b.is_zero():
                A.add_multiple_of_column(i,row - 1, b**t)
                A.add_multiple_of_row(i, row - 1, b)
                D.add_multiple_of_row(i, row - 1, b)

        if row == n - 1:
            break

    # Count how many variables have been used
    if row == n - 1:
        if not A[n - 1, n - 1].is_zero():
            r = n
        else:
            r = n - 1

    # Normalize diagonal elements to 1
    for i in range(r):
        a = A[i, i]
        if not a.is_one():
            # Find an element `b` such that `b*b^t = b^(t+1) = a`
            b = conj_square_root(a)
            D.rescale_row(i, 1 / b)

    return D

@tscrim I've rebased this branch on #39200. It has been updated to call U.cholesky(extended=True).